Robust Process Scheduling (with Python)

Nominal Model

Necessary imports:


In [1]:
%matplotlib inline
from robust_STN import *

We instantiate the nominal model and solve it using the default solver (Gurobi). The data of the nominal model can be changed from within STN.__init__().


In [2]:
stn = STN()
stn.solve()


Constructing nominal model...
Solving...
Parameter OutputFlag unchanged
   Value: 1   Min: 0   Max: 1   Default: 1
Changed value of parameter QCPDual to 1
   Prev: 0   Min: 0   Max: 1   Default: 0
Optimize a model with 1060 rows, 449 columns and 5806 nonzeros
Coefficient statistics:
  Matrix range    [1e-01, 2e+02]
  Objective range [2e-01, 1e+01]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+100]
Found heuristic solution: objective 0
Presolve removed 957 rows and 358 columns
Presolve time: 0.01s
Presolved: 103 rows, 91 columns, 775 nonzeros
Variable types: 51 continuous, 40 integer (40 binary)

Root relaxation: objective -2.908360e+03, 93 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -2908.3600    0   10    0.00000 -2908.3600      -     -    0s
H    0     0                    -2572.000000 -2908.3600  13.1%     -    0s
     0     0 -2889.9539    0    9 -2572.0000 -2889.9539  12.4%     -    0s
     0     0 -2874.5826    0   14 -2572.0000 -2874.5826  11.8%     -    0s
H    0     0                    -2708.000000 -2874.5826  6.15%     -    0s
     0     0 -2871.3632    0   16 -2708.0000 -2871.3632  6.03%     -    0s
     0     0 -2870.2826    0   16 -2708.0000 -2870.2826  5.99%     -    0s
     0     0 -2870.2086    0   17 -2708.0000 -2870.2086  5.99%     -    0s
     0     0 -2870.2086    0   17 -2708.0000 -2870.2086  5.99%     -    0s
     0     2 -2870.2086    0   17 -2708.0000 -2870.2086  5.99%     -    0s
*   18     0               6    -2744.375000 -2866.1550  4.44%   5.8    0s

Cutting planes:
  Gomory: 3
  MIR: 2
  Flow cover: 1

Explored 33 nodes (321 simplex iterations) in 0.06 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective -2.744375000000e+03, best bound -2.744375000000e+03, gap 0.0%
Out[2]:
-2744.3750000000055

Plot the nominal schedule:


In [3]:
stn.plot_schedule()


Attained objective:


In [4]:
stn.model.value


Out[4]:
-2744.3750000000055

Robust Model

A robust STN is instantiated with a reference to a particular nominal STN (e.g., the one created in the previous section).


In [5]:
rSTN = robust_STN(stn)

The module robust_STN.py provides two functions to generate appropriate uncertainty sets:

  • robust_STN.build_uncertainty_set_for_time_delay(units=(0,), tasks=(0,), delay=1, from_t=0, to_t=T)
  • robust_STN.build_uncertainty_set_for_unit_swap(from_unit=2, to_unit=1, tasks=(1,), from_t=0, to_t=T)

Below some examples of how they can be used. These are sufficient to reproduce the results in [1].

Delay Example #1: Heater Delay

We first produce a robust schedule that is immune to possible time delays of the heater (unit=0), when performing heating (task=0), of at most 1 time step (delay=1), which can happen anytime in the scheduling horizon (from_t and to_t assigned default values)


In [6]:
rSTN.W = rSTN.build_uncertainty_set_for_time_delay(units=(0,), tasks=(0,), delay=1)
rSTN.solve()


Constructing robust model...
Solving...
Parameter OutputFlag unchanged
   Value: 1   Min: 0   Max: 1   Default: 1
Changed value of parameter QCPDual to 1
   Prev: 0   Min: 0   Max: 1   Default: 0
Optimize a model with 84422 rows, 69909 columns and 209228 nonzeros
Coefficient statistics:
  Matrix range    [1e-01, 1e+05]
  Objective range [2e-01, 1e+01]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+100]
Found heuristic solution: objective 0
Presolve removed 75604 rows and 64943 columns
Presolve time: 0.44s
Presolved: 8818 rows, 4966 columns, 47405 nonzeros
Variable types: 4895 continuous, 71 integer (71 binary)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -3193.9262    0   18    0.00000 -3193.9262      -     -    0s
H    0     0                    -1965.750000 -3193.9262  62.5%     -    0s
     0     0 -3083.0600    0   25 -1965.7500 -3083.0600  56.8%     -    0s
     0     0 -2908.3600    0   18 -1965.7500 -2908.3600  48.0%     -    1s
     0     0 -2902.1364    0   17 -1965.7500 -2902.1364  47.6%     -    1s
H    0     0                    -2572.000000 -2902.1364  12.8%     -    1s
     0     0 -2901.3522    0   26 -2572.0000 -2901.3522  12.8%     -    1s
H    0     0                    -2708.000000 -2901.3522  7.14%     -    1s
     0     0 -2899.0293    0   24 -2708.0000 -2899.0293  7.05%     -    1s
     0     0 -2898.6051    0   27 -2708.0000 -2898.6051  7.04%     -    1s
     0     0 -2898.5181    0   31 -2708.0000 -2898.5181  7.04%     -    1s
     0     0 -2898.1193    0   24 -2708.0000 -2898.1193  7.02%     -    1s
     0     0 -2897.5287    0   30 -2708.0000 -2897.5287  7.00%     -    1s
     0     0 -2897.5234    0   32 -2708.0000 -2897.5234  7.00%     -    1s
     0     0 -2897.3359    0   28 -2708.0000 -2897.3359  6.99%     -    1s
     0     0 -2897.1707    0   29 -2708.0000 -2897.1707  6.99%     -    1s
     0     0 -2896.8947    0   23 -2708.0000 -2896.8947  6.98%     -    1s
     0     0 -2896.8947    0   25 -2708.0000 -2896.8947  6.98%     -    1s
     0     0 -2896.8748    0   27 -2708.0000 -2896.8748  6.97%     -    2s
     0     0 -2896.8540    0   29 -2708.0000 -2896.8540  6.97%     -    2s
     0     0 -2896.8307    0   25 -2708.0000 -2896.8307  6.97%     -    2s
     0     0 -2896.8307    0   24 -2708.0000 -2896.8307  6.97%     -    2s
     0     2 -2896.8307    0   22 -2708.0000 -2896.8307  6.97%     -    2s
*   45     7              15    -2744.375000 -2847.8872  3.77%  30.0    2s

Cutting planes:
  Gomory: 3
  Implied bound: 34
  Clique: 3
  MIR: 21
  Flow cover: 33
  Zero half: 2

Explored 59 nodes (14103 simplex iterations) in 2.98 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective -2.744375000000e+03, best bound -2.744375000000e+03, gap 0.0%

Plot the robust schedule


In [7]:
rSTN.plot_schedule()


We can then simulate an uncertain event, and see how the schedule has to be adapted to accommodate it.

  • The important feature to observe is that events are accommodated without disrupting the core of the schedule (assignments), but just by regulating processed quantities.
  • Note also that no more heavy computations are required to calculate the adapted schedule.

In [9]:
rSTN.simulate_uncertain_event(event=[0,0,0], column=1)
rSTN.simulate_uncertain_event(event=[0,0,5], column=1)


Events are indexed according to the decision they affect. The column parameter is used to simulate appropriate delays, when they can be longer than 1 time step (column=2 would simulate a delay by two time steps). Check the implementation of robust_STN.simulate_uncertain_event() for more details on how the parameters are to be understood.

Delay Example #2: Reactor 2, delay when processing Reaction 1 or Reaction 3


In [10]:
rSTN.W = rSTN.build_uncertainty_set_for_time_delay(units=(2,), tasks=(1,3), delay=1)
rSTN.solve()


Constructing robust model...
Solving...
Parameter OutputFlag unchanged
   Value: 1   Min: 0   Max: 1   Default: 1
Changed value of parameter QCPDual to 1
   Prev: 0   Min: 0   Max: 1   Default: 0
Optimize a model with 145784 rows, 88989 columns and 390650 nonzeros
Coefficient statistics:
  Matrix range    [1e-01, 1e+05]
  Objective range [2e-01, 1e+01]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+100]
Found heuristic solution: objective 0
Presolve removed 128952 rows and 79523 columns
Presolve time: 1.07s
Presolved: 16832 rows, 9466 columns, 92061 nonzeros
Variable types: 9405 continuous, 61 integer (61 binary)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -2969.1200    0    9    0.00000 -2969.1200      -     -    1s
H    0     0                    -2385.666667 -2969.1200  24.5%     -    1s
     0     0 -2694.3418    0   11 -2385.6667 -2694.3418  12.9%     -    1s
     0     0 -2688.0633    0    9 -2385.6667 -2688.0633  12.7%     -    1s
     0     0 -2688.0633    0   11 -2385.6667 -2688.0633  12.7%     -    1s
     0     0 -2687.9192    0   15 -2385.6667 -2687.9192  12.7%     -    2s
     0     0 -2687.5540    0    9 -2385.6667 -2687.5540  12.7%     -    2s
     0     0 -2686.2644    0   11 -2385.6667 -2686.2644  12.6%     -    2s
     0     0 -2684.7485    0   18 -2385.6667 -2684.7485  12.5%     -    2s
     0     0 -2684.6243    0   19 -2385.6667 -2684.6243  12.5%     -    2s
     0     0 -2684.4342    0   19 -2385.6667 -2684.4342  12.5%     -    2s
     0     0 -2682.8788    0   10 -2385.6667 -2682.8788  12.5%     -    2s
     0     0 -2682.8788    0    9 -2385.6667 -2682.8788  12.5%     -    2s
     0     2 -2682.8788    0    9 -2385.6667 -2682.8788  12.5%     -    2s
*   24     7              16    -2412.000000 -2565.1111  6.35%  82.2    3s
*   49     2               8    -2547.750000 -2553.1111  0.21%  67.2    3s

Cutting planes:
  Gomory: 1
  Implied bound: 15
  MIR: 19
  Flow cover: 11
  Zero half: 2

Explored 62 nodes (5914 simplex iterations) in 3.54 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective -2.547750000000e+03, best bound -2.547750000000e+03, gap 0.0%

In [12]:
rSTN.plot_schedule()
rSTN.simulate_uncertain_event(event=[2,1,0], column=1)
rSTN.simulate_uncertain_event(event=[2,1,5], column=1)


We see that the optimizer does not schedule any reaction 3 on reactor 2. In the above plots we can also inspect how the schedule has to be adjusted when the delay does occur.

Unit Swap Example: swap processing of Reaction 2, from Reactor 1 to Reactor 2, starting from t = 4h


In [13]:
rSTN.W = rSTN.build_uncertainty_set_for_unit_swap(from_unit=1, to_unit=2, tasks=(2,), from_t=4)
rSTN.solve()


Constructing robust model...
Solving...
Parameter OutputFlag unchanged
   Value: 1   Min: 0   Max: 1   Default: 1
Changed value of parameter QCPDual to 1
   Prev: 0   Min: 0   Max: 1   Default: 0
Optimize a model with 63968 rows, 63549 columns and 148754 nonzeros
Coefficient statistics:
  Matrix range    [1e-01, 1e+05]
  Objective range [2e-01, 1e+01]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+100]
Found heuristic solution: objective 0
Presolve removed 62018 rows and 62512 columns
Presolve time: 0.18s
Presolved: 1950 rows, 1037 columns, 9027 nonzeros
Variable types: 988 continuous, 49 integer (49 binary)

Root relaxation: objective -3.142348e+03, 468 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -3142.3480    0   18    0.00000 -3142.3480      -     -    0s
H    0     0                    -2043.000000 -3142.3480  53.8%     -    0s
H    0     0                    -2193.500000 -3142.3480  43.3%     -    0s
     0     0 -2952.6741    0   24 -2193.5000 -2952.6741  34.6%     -    0s
H    0     0                    -2264.416667 -2952.6741  30.4%     -    0s
     0     0 -2748.2120    0   22 -2264.4167 -2748.2120  21.4%     -    0s
     0     0 -2638.1570    0   10 -2264.4167 -2638.1570  16.5%     -    0s
H    0     0                    -2361.416667 -2638.1570  11.7%     -    0s
     0     0 -2627.1500    0   15 -2361.4167 -2627.1500  11.3%     -    0s
     0     0 -2598.9685    0   13 -2361.4167 -2598.9685  10.1%     -    0s
     0     0 -2598.7742    0   14 -2361.4167 -2598.7742  10.1%     -    0s
     0     0 -2598.2411    0   16 -2361.4167 -2598.2411  10.0%     -    0s
     0     0 -2598.0755    0   16 -2361.4167 -2598.0755  10.0%     -    0s
H    0     0                    -2382.635417 -2598.0755  9.04%     -    0s
     0     2 -2598.0755    0   16 -2382.6354 -2598.0755  9.04%     -    0s
*   14     0              10    -2513.750000 -2521.5981  0.31%  32.6    0s

Cutting planes:
  Gomory: 1
  Cover: 2
  Implied bound: 40
  Clique: 4
  MIR: 9
  Flow cover: 19
  GUB cover: 4
  Zero half: 4

Explored 29 nodes (1489 simplex iterations) in 0.55 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective -2.513750000000e+03, best bound -2.513750000000e+03, gap 0.0%

And wen can again check the resulting schedules (remember events are counted from from_t):


In [14]:
rSTN.plot_schedule()
rSTN.simulate_uncertain_event(event=[1,2,4])


Observe how in the revised schedule, the core is maintained, but the batch sizes in the remainder of the schedule are changed (reduced) due to the unit swap to a unit which is smaller.

Notes

1) The code pertaining to the nominal model is within STN.py. It contains the data used in the STN.__init__() method, which now replicates the results from the original paper [2]. It also contains the functions used to construct the 'core' constraints.
2) To run this locally, you need cvxpy and gurobipy (academic lisenses for gurobi can be obtained)
3) Most of the increase in computation time from the nominal model is due to the computation of the recourse matrix. If you do not need recourse, or your problem doesn't have/allow the adjustment of the continuous variables when uncertain outcomes occur, you can use the RC provided in Corollary 2.4 in [1]; the resulting robust cuonterpart then is exactly as expensive as the nominal model.
4) The reported solve times on this notebook may be different from those in the paper; they were solved on a different machine (the results in thie notebook are based on the performance on a laptop).

References

  1. R. Vujanic, P. Goulart, M. Morari, Robust Optimization of Schedules Affected by Uncertain Events, July 2015, submitted.
  2. E. Kondili, C.C. Pantelides, and R.W.H. Sargent, A general algorithm for short-term scheduling of batch operations - I. MILP formulation., Computers & Chemical Engineering 17 (1993), no. 2, 211–227.